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Abstract 

A  three  dimensional  model  of  a  proton  exchange  membrane  fuel  cell  (PEMFC)  operating  with  a  polybenzimidazole  (PBI)  membrane  is  presented. 
This  model  is  an  improvement  on  the  previous  three  models  developed  for  this  type  of  PEM  fuel  cell.  The  model  accounts  for  all  transport  and 
polarization  phenomena,  and  the  results  compare  well  with  published  experimental  data  for  the  same  sets  of  operating  conditions.  The  model  predicts 
the  oxygen  depletion,  which  occurs  in  the  catalyst  area  under  the  ribs,  and  a  temperature  rise  of  up  to  20  K  at  a  power  density  of  1000  W  m-2  can 
be  expected  depending  on  operating  conditions. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Recently  there  has  been  a  growing  interest  in  high  tempera¬ 
ture  (120-150  °C),  low  humidity  (25-50%  RH)  proton  exchange 
membrane  fuel  cells  (PEMFCs).  The  US  Department  of  Energy 
(DOE)  has  planned  to  invest  US$  17.5  million  into  research  in 
such  fuel  cells  over  the  next  5  years  [1]. 

Typical  PEMFCs,  equipped  with  a  Nafion®  membrane,  oper¬ 
ate  at  temperatures  below  the  boiling  point  of  water  (50-80  °C), 
since  Nafion®  requires  a  high  water  content  to  effectively  con¬ 
duct  protons.  To  maintain  these  low  temperatures  under  vehicle 
driving  conditions,  especially  at  peak  power,  requires  over-sized 
cooling  equipment.  External  humidification  is  also  required 
to  prevent  dehydration  of  the  Nafion®  membrane,  and  this 
adds  volume,  weight  and  complexity  to  the  fuel  cell  system 
[1]. 

At  higher  temperatures,  water  exists  primarily  in  the  vapor 
phase,  so  transport  limitations  associated  with  the  presence  of 
liquid  water  are  precluded.  In  addition,  carbon  monoxide  (CO) 
adsorption  onto  the  platinum  (Pt)  catalyst  particles  becomes  less 
pronounced  at  elevated  temperatures.  Thus,  its  poisonous  effects 
on  fuel  cell  performance  decrease  with  temperature.  Also  with 
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high  temperature  membranes,  which  are  not  dependant  on  water 
for  proton  conductivity,  external  humidification  is  not  required. 
This  reduces  overall  system  cost,  weight  and  size. 

Proton  exchange  membranes  fall  into  three  general  cate¬ 
gories:  (1)  per-fluorinated  ionomer  membranes;  (2)  partially 
per-fluorinated  membranes;  (3)  non-per-fluorinated  membranes 
[2].  Nafion®  is  an  example  of  a  per-fluorinated  ionomer  mem¬ 
brane,  whereas  polybenzimidazole  (PBI)  is  an  example  of  a 
non-per-fluorinated  membrane.  PBI  was  first  suggested  for  use 
in  fuel  cells  by  Wainright  et  al.  [3].  PBI  has  good  mechanical 
strength,  high  chemical  and  thermal  stability  at  high  temper¬ 
atures,  and  it  attains  its  ionic  conductivity  when  doped  with  a 
strong  acid  such  as  phosphoric  or  sulfuric  acid.  PBI  is  also  highly 
impermeable  to  methanol  at  high  temperatures,  and  has  a  water 
drag  coefficient  of  nearly  zero,  which  alleviates  cathode  flooding 
and  membrane  dehydration. 

The  ionic  conductivity  of  PBI  varies  significantly  depending 
on  temperature,  the  level  of  acid  doping,  method  of  preparation 
and  composite  materials  used,  if  any.  Typical  conductivity  values 
reported  for  PBI  doped  with  phosphoric  acid  range  from  0.2  to 
6.8  Sm-1  [3-7];  5-6  S  m-1  when  doped  with  sulfuric  acid  [8,9]; 
and  9  S  m_1  when  doped  with  potassium  hydroxide  [7].  A  com¬ 
posite  of  PBI  and  ZrP  has  exhibited  a  conductivity  of  9.6  S  m-1 
[4],  while  the  highest  conductivity  value  reported  in  the  litera¬ 
ture  is  13  S  m_1  at  160  °C  for  PBI  doped  with  1300-1600  mol% 
phosphoric  acid  [7] .  There  is  still  a  debate  over  the  dependence 
of  PBI  conductivity  on  humidity,  whether  any  such  dependence 
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Nomenclature 

a  effective  surface  area  (m-  1 ) 

cp  specific  heat  capacity  at  constant  pressure 
(Jkg-1  K-1) 

Dij  diffusivity  of  gas  pair  i-j  (m2  s-1) 

Dj  effective  thermal  diffusivity  (m2  s-1) 

F  Faraday  constant,  96,487  C  mol-1 

i  current  density  (A  m-2) 

i{)  exchange  current  density  (A  m-2) 

j  reaction  rate  (A  m-  3 ) 

k  thermal  conductivity  (W  m-1  K-1) 

kp  gas  permeability  (m2) 

m  mass  fraction 

M  molar  mass  (kg  mol- 1 ) 

P  pressure  (Pa) 

R  universal  gas  constant,  8.3143  J  mol-1  K-1 
S  source,  entropy 

T  Temperature  (K) 

u  velocity  (ms-1) 

V  potential  (V) 

v  mole  fraction 

Greek  letters 

a  thermal  diffusivity  (m2  s  - 1 ) ,  charge  transfer  coef¬ 

ficient 

y  concentration  parameter 

s  porosity 

/j  dynamic  viscosity  (Pa  s) 

p  density  (kg  m-3) 

a  electrical  or  ionic  conductivity  (S  m-1) 

0  electrical  or  ionic  potential  (V) 

Subscripts 

e  electrolyte  phase 

f  fluid  property 

i,  j  species  i,  j 

m  mass  quantity 

s  solid  phase 

T  thermal 

Superscripts 
d  dispersion 

eff  effective 


exists,  and  if  so  the  nature  of  it  [10-12].  PBI  has  been  success¬ 
fully  tested  in  PEMFCs  using  hydrogen  [8,9,13,14],  methanol 
[15]  and  other  alcohols  [16,17].  Maximum  power  densities  are 
in  the  order  of  1000  W  m-2  using  methanol  as  the  fuel  [15]  and 
4000-9000  W  m-2  using  hydrogen  as  the  fuel  [8,9]. 

Most  of  the  work  done  in  fuel  cell  modeling  over  the  past 
15  years  has  been  for  Nation®  membranes  [18,19].  The  only 
attempts  thus  far  to  model  fuel  cells  equipped  with  PBI  mem¬ 
branes  have  been  conducted  by  the  present  authors  [20-22] .  This 
paper  builds  on  our  previous  three  models,  where  we  studied  ID 
and  2D  effects.  This  paper  presents  a  3D  non-isothermal  model 


of  a  PBI  fuel  cell.  The  model  describes  all  transport  and  polar¬ 
ization  phenomena,  and  accounts  for  rib  effects  and  the  effect  of 
gas  channel  flow  rates.  Numerical  results  for  a  hydrogen/oxygen 
and  a  hydrogen/air  fuel  cell  are  compared  to  the  experimental 
results  published  by  Wang  et  al.  [13]. 

2.  Model  development 

2.7.  Computational  domain 

Fig.  1  shows  the  computational  domain  used  in  this  model, 
as  well  as  the  grid  network.  In  the  x-direction,  the  domain  spans 
the  anode  gas  channel  and  rib  to  the  cathode  gas  channel  and 
rib.  It  includes  the  membrane  electrode  assembly  (ME A),  with 
catalyst  layers  treated  as  finite  sized  regions,  as  opposed  to  van¬ 
ishing  interfaces.  In  the  y-direction,  the  domain  spans  half  of 
one  channel  to  half  of  an  adjacent  rib.  In  the  work  by  Wang  et  al. 
[13],  their  fuel  cell  testing  apparatus  used  five  straight  gas  chan¬ 
nels.  To  simplify  the  computer  memory  requirements,  analysis 
was  only  performed  for  one  channel,  symmetry  being  assumed 
at  the  mid-points  of  each  channel  and  rib.  In  the  z-direction,  the 
domain  spans  the  entire  length  of  the  gas  channel  passing  over 
the  active  MEA,  which  is  1  cm  x  1  cm  in  cross  section.  Thus, 
the  length  of  the  channels  in  the  domain  is  1  cm.  This  dimen¬ 
sion  is  an  order  of  magnitude  larger  than  those  in  the  x-  and 
y-directions.  This  geometric  aspect  ratio  is  not  obvious  in  Fig.  1 
for  the  purpose  of  making  the  figure  clearer. 

2.2.  Assumptions 

For  high  temperature  fuel  cells,  water  is  expected  to  exist  only 
in  the  vapor  phase;  therefore,  single  phase  flow  is  assumed.  Indi¬ 
vidual  gases  as  well  as  the  gas  mixtures  are  assumed  to  behave 
ideally.  The  porous  media  and  catalyst  regions  are  assumed  to 
be  isotropic  and  macro-homogeneous,  and  the  PBI  membrane 
is  assumed  to  be  impermeable  to  gas  flow.  Also  because  of  the 
typically  low  fluid  velocities,  the  flow  is  treated  as  laminar. 

The  electrochemical  reactions  are  assumed  to  be  gas  phase 
reactions.  The  present  model  does  not  take  into  account  reac¬ 
tant  gases  dissolving  in  the  electrolyte  at  the  catalyst  layer.  This 
would  require  two-phase  modeling,  which  is  the  subject  of  future 


Anode  Gas  Channel  Membrane  Cathode  Gas  Channel 


Fig.  1.  3D  computational  domain  and  finite  element  grid. 
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work.  It  is  expected  that  a  two-phase  model,  which  takes  into 
account  gas  solubility  in  the  catalyst  layers,  would  be  better 
able  to  predict  mass  transport  limitations.  This  model  also  does 
not  take  into  account  the  effect  of  feed  gas  humidification  on 
membrane  conductivity  and  electro-kinetic  activity,  since  the 
mechanisms  governing  these  phenomena  are  not  as  yet  fully 
understood. 

2.3.  Governing  equations 

V  .{pu)  =  Sm  =  YJSi  (1) 


The  charge  conservation  equations  in  the  electrolyte  and  solid 


phases  are  given  by, 

V  •  is  =  V .  (~afy<ps)  =  -j 

(8) 

v  •  ie  =  v  .  (-a eeffV0e)  =  +j 

(9) 

•  -ref 

J  =  ai0 


XjP 
P ref 


aF 


—  exp 


(10) 


The  continuity  equation  (1)  has  a  non-zero  term  on  the  RHS, 
which  represents  a  solid-fluid  phase  change  in  the  catalyst  lay¬ 
ers.  For  a  PBI  membrane,  H+  ions  are  transported  across  the 
membrane  in  a  solid  state  via  a  Grothus  mechanism  [10-12].  So 
there  is  a  loss  of  fluid  mass  at  the  anode  catalyst  layer  and  a  gain 
in  fluid  mass  at  the  cathode  catalyst  layer,  as  shown  in  the  half 
cell  equations  (2)  and  (3). 

H2  (g)  — >  2H+  (s)  +  2e_  (anode)  (2) 

5O2  (g)  +  2H+  (s)  +  2e_  ->  H2O  (g)  (cathode)  (3) 

Total  mass,  however,  is  always  conserved  at  every  point;  and 
overall  fluid  mass  is  conserved  since  the  H+,  formed  at  the  anode, 
reconverts  to  fluid  mass  at  the  cathode.  But  at  the  catalyst  lay¬ 
ers,  the  fluid  flow  is  non-conservative.  Across  the  membrane, 
the  fluid  velocity  is  zero  since  the  PBI  membrane  is  virtually 
impermeable  to  gas  flow. 

The  Navier-Stokes  equations  govern  momentum  transfer  in 
the  gas  channels  and  the  electrodes. 

pu  •  Vw  =  -VP  +  V  •  (gVw)  —  —u  (4) 

kv 

For  the  porous  diffuser  and  catalyst  regions,  the  source  term 
corresponding  to  Darcy’s  law  becomes  dominant  over  the  inertia 
and  viscous  terms.  This  source  term  is  non-existent  in  the  gas 
channels,  where  the  flow  is  primarily  inertial. 

Considering  that  a  non-conservative  form  of  the  continuity 
equation  is  used,  the  non-conservative  Stefan-Maxwell  equa¬ 
tions  must  be  used  for  gas  species  conservation,  which  were 
previously  derived  [21,22]. 


pu  •  Vm;  =  V  • 


N 


psmi 


fvm,  +  nii 


T  Si  int  3Vn 


The  last  term  on  the  RHS  accounts  for  the  non-conservative 
nature  of  the  fluid  flow. 

The  energy  equation  accounts  for  convection,  conduction  and 
heat  generation  due  to  ohmic  or  Joule  heating,  as  well  as  heat  of 
reaction. 


pcpU  •  VP  =  pcpW  •  (DtVT)  +  St 

(6) 

St  =  ^ohm  T  SVXn 

(7) 

The  exchange  current  density,  io,  is  defined  in  terms  of  the 
active  catalyst  surface  area.  The  effective  surface  area,  a ,  is 
defined  as  the  ratio  of  the  total  active  catalyst  surface  area  to  the 
total  catalyst  region  volume,  and  thus  takes  into  account  surface 
roughness  in  the  catalyst  layer.  The  Butler- Volmer  equation  (10) 
states  that  the  rate  of  electrochemical  reaction  is  driven  by  the 
difference  in  potential  between  the  two  respective  phases,  and 
is  affected  by  the  concentration  of  reactants  at  the  catalyst  sites. 
The  local  reaction  rate,  j,  defines  the  rate  of  transfer  of  solid  state 
current  to  electrolyte  phase  current  in  the  anode  catalyst  layer, 
and  vice  versa  in  the  cathode  catalyst  layer. 


2.4.  Source  terms 


The  source  terms  depend  on  the  reaction  rate  in  each  catalyst 
layer  and  the  stoichiometry  of  the  overall  cell  reaction.  With 
the  exception  of  the  ohmic  heating  term,  all  other  source  terms 
are  only  non-zero  in  the  catalyst  layers.  The  reaction  rate,  j,  is 
positive  in  the  anode  catalyst  layer  and  negative  in  the  cathode. 
Thus,  the  following  equations  define  the  local  rates  of  oxygen 
consumption,  water  vapor  production,  hydrogen  consumption, 
reactive  heat  generation  and  ohmic  heat  generation: 


So2  =  j 


Sr2  =  —  j 


3ohm  — 


^°2  (kg  m  3  s  1 ) 

AF  5 

(11) 

.Mh20  ,,  _3  _K 

J  OZ7  (kgm  s  ) 

2  F 

(12) 

~j  2f  (kgm  3s  *) 

(13) 

-j(fc-<P s-^r)  (Wm“3) 

(14) 

CTeSff  +  aeeff(Wm_3) 

(15) 

2.5.  Constitutive  relations 

This  section  provides  a  list  of  constitutive  relations  necessary 
to  correct  the  plain  media  fluid  and  material  properties  for  the 
porous  media.  It  also  gives  the  ideal  gas  relations  used  in  this 
work.  The  geometric  average  of  the  solid  and  fluid  thermal  con¬ 
ductivity  is  used  to  determine  the  effective  thermal  conductivity 
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[23]. 

kefi  =  As\1_£ 

kt  \kf) 


(16) 


The  thermal  dispersion  diffusivity  is  given  by  the  following  rela¬ 
tionship,  valid  for  low  Peclet  number  flows  [24]: 


Oif 


\(k8/kf)  +  2j 


(17) 


The  total  effective  thermal  diffusivity  is  given  by  Kaviany  [25], 


Dj  kcft  sDd 

=  T  1 

off  kf  off 


(18) 


This  value  depends  on  solid  and  fluid  properties  and  composi¬ 
tions,  as  well  as  porous  media  characteristics.  The  effective  gas 
diffusivity  is  related  to  the  plain  media  diffusivity  [26], 


Fj  2 £ 

Aj  3  —  s 


(19) 


The  electrical  and  ionic  conductivities  are  corrected  for  porous 
media  using  the  Bruggemann  correlation  [27],  where  £phase  is 
the  volume  fraction  occupied  by  the  phase  through  which  the 
respective  current  flows. 

_eff _  1.5 

°  (2d) 


The  gas  mixture  (mass)  density  and  molar  mass  are  derived  from 
the  ideal  gas  law. 


PM 

~RT 


(21) 


M  = 


(22) 


All  other  fluid  properties  are  mass  averaged  for  the  gas  mix¬ 
ture.  The  plain  media  gas  pair  diffusivity  is  assumed  constant  at 
a  given  temperature  and  pressure,  and  independent  of  composi¬ 
tion. 


2.(5.  Boundary  conditions 

The  unified  (single  domain)  approach  is  used  in  this  model. 
Thus,  boundary  conditions  are  required  at  the  ends  of  the 
domain.  In  the  experiments  by  Wang  et  al.  [13],  the  feed  gases 
are  humidified  at  28  °C  and  supplied  to  the  fuel  cell  at  150  °C 
and  atmospheric  pressure.  Their  PBI  membrane  after  doping 
with  5  M  phosphoric  acid  had  a  thickness  80  jim.  They  used 
a  commercial  E-TEK  electrode  (0.5  mg  cm-2  Pt)  at  the  anode, 
and  a  home  made  cathode  with  2  mg  cm-2  Pt  loading.  They 
did  not,  however,  report  the  supply  gas  flow  rates  used  in  their 
experiments.  So  in  the  computations,  typical  values  were  used. 

In  the  extremes  of  the  domain  in  the  y-direction,  Neumann 
boundary  conditions  are  specified,  since  symmetry  conditions 
were  assumed  in  the  selection  of  the  domain.  In  the  extremes  in 
the  ^-direction,  the  flow  inlet  and  outlet  conditions  are  specified. 
The  overall  flow  rate  is  specified  at  the  inlets  of  each  channel, 
while  atmospheric  pressure  is  specified  at  the  outlet.  For  the 


species  and  thermal  equations,  the  concentration  and  temper¬ 
ature,  respectively,  are  specified  at  the  inlet,  while  convective 
flux  boundary  conditions  are  specified  at  the  outlets.  In  the  x- 
direction,  insulation  conditions  are  specified  for  the  mass  flow 
and  species  flow  at  the  membrane/catalyst  interfaces  since  the 
membrane  is  impermeable  to  gas  flow.  The  solid  phase  poten¬ 
tials  are  specified  at  the  ribs  in  the  anode  and  cathode.  Zero  is 
chosen  as  a  reference  at  the  anode,  while  Vrev  —  VCeii  (the  differ¬ 
ence  between  the  reversible  cell  potential  and  the  operating  cell 
potential)  is  specified  at  the  cathode.  This  quantity  represents 
the  total  cell  polarization.  Thus,  in  this  model,  the  cell  potential 
is  stated  and  the  cell  current  density  computed. 

3.  Results  and  discussion 

The  governing  equations  are  solved  using  FEMFAB®  3. li, 
utilizing  a  finite  element  method  to  solve  the  system  of  coupled 
partial  differential  equations.  The  computations  were  performed 
on  a  Windows  PC  with  a  1  GB,  3.4  GHz,  32  bit,  Pentium  4  pro¬ 
cessor.  The  finite  element  method  entailed  3190  second  order 
tetrahedral  elements,  which  for  10  independent  variables,  consti¬ 
tuted  a  25,138  degree  of  freedom  problem  with  a  computational 
error  O (h3).  The  computational  grid  is  shown  in  Fig.  1.  Com¬ 
putation  times  were  up  to  5  min  for  a  single  run.  A  run  consists 
of  computing  the  entire  solution  for  one  specified  cell  voltage. 
Once  this  initial  solution  is  obtained,  it  can  be  used  as  an  initial 
guess  to  compute  the  solutions  for  other  cell  voltages.  Subse¬ 
quent  runs  take  approximately  3  min  to  solve.  Thus,  all  the  data 
for  an  IV  curve  can  be  obtained  in  approximately  30  min. 

This  size  of  grid  was  considered  the  best  compromise 
between  solution  time  and  accuracy.  Fig.  2  shows  the  tests  for 
grid  independence.  The  variable  which  showed  the  most  signif¬ 
icant  variation  with  grid  size  was  found  to  be  the  cell  current 
density,  especially  at  limiting  conditions.  So  this  quantity  was 
shown  in  Fig.  2  versus  the  number  of  grid  elements.  Taking 
into  consideration  the  scale  used  in  the  figure,  even  coarse  grids 
allowed  for  an  error  of  less  than  0.5%.  However,  the  solution 
time  increased  significantly  when  the  grid  was  refined.  A  grid 
consisting  of  3190  elements  required  approximately  5  min  of 
solution  time  and  gave  an  error  of  approximately  0.2%,  while 
a  grid  of  8600  elements  required  over  20  min  of  solution  time 


Fig.  2.  Grid  independence  test. 
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Table  1 


Membrane  and  diffuser  properties 


Membrane  [13] 

Thermal  conductivity,  km  (W  m-1  K_1)  [31] 

40 

Ionic  conductivity,  crm  (S  m-1) 

1.87 

Thickness  (m) 

8  x  10“5 

Cross  section  (m2) 

0.01  X  0.01 

Diffuser  (graphite)  [28] 

Thermal  conductivity,  kd  (W  m-1  K-1)  [32] 

1.15 

Electrical  conductivity,  crd  (S  m-1) 

120 

Porosity,  e 

0.4 

Permeability,  kv  (m2) 

1.8  x  1CT11 

Thickness  (m) 

2.6  x  10“4 

Channels  [13] 

Channel  length  (m) 

0.01 

Channel  cross  section  (deduced)  (m2) 

0.001  X  0.001 

Rib  cross  section  (deduced)  (m2) 

0.001  X  0.001 

Number  of  channels 

5 

Table  2 

Catalyst  layer  properties  at  423  K 

Anode  Cathode 

a(*f(Am-3)  [30] 

108  0.02 

a  [28] 

0.5  2 

V  [28] 

0.5  1 

Aef  (P&) 

1.013  x  105 

s  (estimated) 

0.2 

Fraction  of  membrane  phase  in  catalyst  layer,  £mx  [28] 

0.4 

Catalyst  layer  thickness  (m) 

5  x  10“5 

Table  3 

Fluid  properties  at  423  K  [33] 

Oxygen  Nitrogen 

Water  vapor  Hydrogen 

M  (kg  mol  1 ) 

32  x 10“3 

28.16  x  10“3 

18  x 10“3 

2  x  10“3 

k  (W  m-1  K-1) 

0.036 

0.034 

0.030 

0.239 

Cp  (Jkg-1  K-1) 

956 

1050 

1980 

14500 

U  (Pa  s) 

27  x 10“6 

23  x  10“6 

14  x 10“6 

1.1  x  10“ 

a  (m2  s-1) 

44.4  x  10“6 

41.0  x  10“6 

30.8  x  10“6 

217  x  10“ 

with  only  a  slight  improvement  in  accuracy.  It  was  concluded 
that  grid  independence  was  attained  with  3190  elements. 

Tables  1-4  give  a  list  of  numerical  values  used  in  the  compu¬ 
tations.  The  gas  diffusers  used  by  Wang  et  al.  [13]  are  the  same 
as  those  used  for  PEMFCs  using  Nation®,  so  the  gas  diffuser 
parameters  are  typical  values  taken  from  the  literature  [28].  In 
the  work  of  Wang  et  al.  [13],  their  fuel  cell  used  a  MEA  of  cross 
section  1  cm2,  as  well  as  flow  field  containing  five  straight  gas 
channels.  The  dimensions  of  the  gas  channels  and  ribs  are  calcu¬ 
lable  4 

Plain  media  gas  pair  diffusivities  at  423  K  and  1  atm  [34] 

Gas  pair  Diffusivity  (m2  s-1) 


o2-h2o 

41.9  x  10“6 

o2-n2 

34.2  X  10“6 

n2-h2o 

49.2  x  10“6 

h2-h2o 

144  x  10“6 

lated  assuming  that  the  channels  are  equally  spaced  and  that  rib 
and  channel  thicknesses  are  approximately  the  same.  The  elec¬ 
trode  kinetics  data  for  phosphoric  acid  doped  PBI  are  reported 
to  be  approximately  the  same  as  that  of  pure  phosphoric  acid 
[29],  so  electro-kinetic  data  for  a  phosphoric  acid  fuel  cell  are 
obtained  from  Choudhury  et  al.  [30]. 

3.1.  Polarization  curves 

Fig.  3  shows  the  polarization  curves  for  both  hydro¬ 
gen/oxygen  and  hydrogen/air  fuel  cells.  The  experimental  data 
are  only  given  by  Wang  et  al.  [13]  for  up  to  2500  A  m-2,  which 
covers  the  activation  and  some  of  the  ohmic  regions  of  the  IV 
curve.  It  is  difficult  to  say,  from  their  data,  what  limiting  current 
density  was  observed.  However,  for  the  ohmic  and  activation 
regions,  the  model  curves  match  the  experimental  polarization 
data  well.  The  model  results  for  air  match  extremely  well  with 
the  experimental  results.  For  oxygen,  there  is  a  good  match 
except  for  the  last  two  data  points,  where  the  models  underesti¬ 
mate  the  performance.  The  last  two  points  in  the  experimental 
data  for  oxygen  are  indeed  strange,  for  it  appears  as  though  there 
is  an  increase  in  membrane  conductivity  as  the  current  density 
increases.  Barring  the  possibility  of  this  being  experimental  error 
(not  allowing  the  cell  to  fully  reach  steady  state  conditions),  it 
is  possible  that  the  higher  production  of  water  vapor  at  higher 
current  densities  results  in  increased  membrane  conductivity,  or 
maybe  even  enhanced  transport  characteristics  in  the  catalyst 
layer  leading  to  more  rapid  electrode  kinetics. 

Fig.  3  also  shows  the  effects  of  varying  the  inlet  gas  flow  rates 
on  polarization  behavior.  For  the  oxygen  cell,  there  is  little  differ¬ 
ence  as  the  flow  rate  increases  from  0.02  to  0. 1  F  min-1 .  A  flow 
rate  of  0.02  F  min-1  is  not  typically  used  in  laboratory  fuel  cells; 
nevertheless,  it  is  included  in  the  computation  solely  for  the  pur¬ 
pose  of  comparison.  For  the  air  cell,  there  is  little  difference  in 
the  polarization  behavior  as  the  flow  rate  varies  between  0. 1  and 
1  Emin-1;  however,  the  variation  is  much  larger  between  0.02 
and  0.1  Emin-1.  The  limiting  current  density  over  this  range 
varies  from  approximately  4500  to  6300  Am-2.  However,  the 
differences  between  these  two  IV  curves  are  primarily  evident  at 
limiting  conditions.  This  suggests  that  the  air  stream  is  unable  to 
supply  oxygen  at  the  required  rate  for  the  electrochemical  reac¬ 
tion  to  take  place.  It  can  be  shown  that  the  rate  of  consumption 


Fig.  3.  Current-voltage  ( IV)  curves — model  predictions  and  experimental  data. 
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Current  Density  (A  nr2) 


Fig.  4.  Comparison  of  TV  curves  for  3D  and  2D  [Cheddie]  models,  flow 
rate  =  3  L  min~ 1 . 


have  a  distinctly  steeper  gradient,  which  implied  a  larger  ohmic 
resistance.  The  reason  for  this  is  that  the  3D  model  accounts 
for  rib  effects.  The  ribs  provide  an  additional  contribution  to  the 
ohmic  overpotential.  But  also  because  of  the  presence  of  ribs, 
there  is  a  greater  mean  electron  path,  meaning  that  the  current 
must  flow  through  a  slightly  greater  distance  to  reach  the  ribs, 
than  if  the  ribs  were  not  present.  This  results  in  a  further  ohmic 
potential  drop.  So  the  combination  of  both  effects  results  in  a 
greater  overall  cell  resistance.  The  consequence  is  a  smaller  lim¬ 
iting  current  density  predicted  by  the  3D  models.  The  differences 
in  limiting  current  density  predicted  by  2D  and  3D  models  are 
between  1000  and  2000  A  m-2.  The  presence  of  ribs  also  results 
in  areas  of  oxygen  depletion,  which  results  in  slower  electrode 
kinetics.  However,  this  did  not  appear  to  have  a  significant  effect 
on  the  TV  curves.  The  major  difference  between  the  2D  and  3D 
curves  appears  to  be  ohmic  in  nature. 


of  oxygen  in  the  electrochemical  reaction  is  given  by  Eq.  (23), 
from  which  at  4500  Am-2,  the  rate  of  oxygen  consumption  is 
0.04  mg  s- 1 .  For  an  inlet  air  flow  rate  of  0.02  L  min- 1 ,  the  mass 
supply  rate  of  oxygen  is  0.06  mg  s-1 . 

„  _  Ae  hAmea^o2 

K  oxygen  — 

Thus,  the  air  stream  is  barely  able  to  supply  the  necessary 
oxygen  requirement,  which  results  in  the  distinct  concentration 
overpotential  region  for  this  flow  rate.  But,  of  note  is  that  this  is 
the  only  case  where  a  distinct  concentration  overpotential  region 
is  seen.  The  reason  there  was  no  similar  concentration  region 
in  the  other  curves  is  because  there  is  no  liquid  water  at  these 
high  temperatures  to  restrict  the  flow  of  reactant  gases  to  the 
catalyst  layer.  There  may  be  other  factors  which  result  in  mass 
limiting  conditions,  such  as  poor  reactant  transport  and  kinetics 
in  the  catalyst  layer.  However,  this  model  does  not  account  for 
such  phenomena,  neither  does  the  experimental  data  provide  any 
indication  of  what  are  the  limiting  current  densities. 

Fig.  4  gives  a  comparison  of  the  TV  curve  for  the  present  3D 
model  and  our  previous  2D  model  [21]  for  a  common  flow  rate  of 
3  F  min-1 .  For  the  3D  model,  the  ohmic  portions  of  the  TV  curve 


3.2.  Transport  characteristics 

Oxygen  depletion  is  a  problem,  which  occurs  when  air  is 
used  as  the  oxidant  because  of  its  lower  oxygen  content  than 
that  of  pure  oxygen.  Figs.  5  and  6  show  the  oxygen  mass  frac¬ 
tion  variations  in  the  cathode  region  for  two  flow  rates,  1.0  and 
0. 1  F  min-1 ,  respectively.  These  figures  apply  at  limiting  condi¬ 
tions  (zero  cell  voltage),  and  for  humidified  air  as  the  oxidant.  It 
shows  that  concentration  effects  are  more  pronounced  at  lower 
flow  rates.  For  both  flow  rates,  there  is  a  region  of  low  oxy¬ 
gen  concentration  in  the  catalyst  layer  under  the  ribs.  There  is 
a  general  decrease  in  concentration  in  the  negative  v-direction, 
since  oxygen  must  diffuse  against  the  bulk  flow.  There  is  also 
a  concentration  decrease  in  the  negative  y-direction  due  to  the 
presence  of  the  ribs.  Oxygen  must  travel  a  greater  distance  to 
diffuse  from  the  channels  to  the  area  under  the  ribs.  This  oxy¬ 
gen  depletion  is  a  slightly  more  pronounced  at  the  lower  flow 
rate,  and  results  in  slow  reaction  kinetics  and  ineffective  cata¬ 
lyst  utilization.  However,  the  biggest  difference  between  the  two 
flow  rates  occurs  in  the  z-direction.  There  is  a  greater  variation 
in  oxygen  concentration  along  the  gas  channel  in  the  direction 
of  flow  (positive  z-direction).  For  the  lower  flow  rate,  a  larger 


OXYGEN  MASS  FRACTION  Max:  0.21 9 

0.2 

0.18 

0.16 

0.14 

0.12 
0.1 

0.08 
Mint:  0.0724 

Fig.  5.  Oxygen  mass  fraction  for  the  humidified  hydrogen/air  cell,  flow  rate  =  1  Lmin-1 . 
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OXYGEN  MASS  FRACTION  Max:  0.21 8 

0.2 

0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 
Min:  0.0547 

Fig.  6.  Oxygen  mass  fraction  for  the  humidified  hydrogen/air  cell,  flow  rate  =  0. 1  L  min-1 . 


TEMPERATURE 


Max:  429.51 1 


Min:  423.105 

Fig.  7.  Temperature  (K)  variation  for  the  humidified  hydrogen/air  cell,  flow  rate  =  1  L  min-1 . 


fraction  of  oxygen  is  consumed  in  the  electrochemical  reactions. 
The  result  is  a  visible  decrease  in  oxygen  mass  fraction  along 
the  channel  at  the  lower  flow  rate  (Fig.  6).  This  variation  is  not 
as  visible  at  the  higher  flow  rate  (Fig.  5). 


Figs.  7  and  8  show  the  temperature  variations  in  the  cell 
for  the  same  two  flow  rates  mentioned  above.  These  pertain 
to  the  air  cell  operating  at  a  cell  voltage  of  0.4  V,  which  is 
close  to  optimum  operation.  They  show  a  temperature  varia- 


TEMPERATURE  Max:  442.024 

442 

440 

438 

436 

434 

432 

430 

428 

Min:  426.193 

Fig.  8.  Temperature  (K)  variation  for  the  humidified  hydrogen/air  cell,  flow  rate  =  0. 1  L  min- 1 . 
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Fig.  9.  Temperature  variation  and  power  density  vs.  cell  voltage. 


tion  in  the  cell  of  7  and  19  K,  respectively,  the  smaller  variation 
for  the  higher  flow  rate.  From  a  simple  mass  balance  calcu¬ 
lation,  it  can  be  shown  that  such  temperature  variations  are 
expected.  The  volumetric  heat  capacity  in  both  channels  is 
approximately  1000Jm-3K-1.  The  total  heat  generation  for 
the  conditions  represented  in  Figs.  7  and  8  is  0.015  W  for  the 
computational  domain  (0.077  W  for  the  entire  fuel  cell  unit, 
recalling  that  there  are  five  gas  channels).  This  heat  genera¬ 
tion  was  obtained  from  the  software  by  integrating  the  local 
heat  generations  over  the  entire  computational  domain.  For  flow 
rates  of  1 .0  and  0. 1  L  min- 1 ,  the  expected  temperature  “pick  up” 
by  the  gas  streams  are  expected  to  be  approximately  2  and  20  K, 
respectively.  This  is  consistent  with  the  temperatures  in  the  gas 
channels  shown  in  Figs.  7  and  8. 

Fig.  7  (and  to  a  less  visible  extent,  Fig.  8)  shows  that  there 
are  temperature  variations  in  all  three  co-ordinate  directions.  The 
z-direction  variation  is  the  result  of  heat  being  absorbed  by  the 
gas  streams.  The  heat  produced  in  the  fuel  cell  is  transported,  via 
convection,  by  the  flow  gases.  The  v-direction  variation  is  due 
to  the  low  ionic  conductivity  of  the  PBI  membrane  compared  to 
the  gas  diffusion  electrodes.  This  results  in  higher  ohmic  heating 
in  the  membrane.  There  is  also  heat  produced,  which  is  assumed 
to  be  greatest  in  the  cathode  catalyst  layer.  This  heat  must  flow 
out  of  the  MEA  in  the  v-direction,  resulting  in  a  temperature 
variation  in  this  direction,  with  the  hot  points  occurring  close  to 
the  membrane.  The  most  significant  y-direction  variation  is  seen 
where  the  gas  channels  and  the  ribs  meet.  There  is  ohmic  heating 
in  the  ribs,  and  negligible  heat  generation  in  the  channels.  As  a 
result,  the  ribs  attain  a  higher  temperature  than  the  channel  bulk 
flow. 

These  results  were  shown  for  a  cell  voltage  of  0.4  V.  Fig.  9 
shows  the  variation  of  maximum  temperature  variation  versus 
cell  voltage  for  the  two  flow  rates  of  interest.  The  figure  also 
shows  the  power  density  variation,  which  peaks  between  0.3  and 
0.4  V.  The  power  density  is  not  significantly  different  for  the  two 
flow  rates,  so  only  one  curve  is  shown.  The  temperature  vari¬ 
ation  increases  exponentially  as  limiting  conditions  approach 
(cell  voltage  tends  to  zero),  where  temperature  increases  of 
20-70  K  are  expected.  The  reason  for  this  is  that  as  the  cell  volt¬ 
age  decreases,  the  current  density  increases  resulting  in  greater 


Fig.  10.  Ohmic  and  reactive  heat  generation  vs.  cell  voltage. 

ohmic  heating.  The  heat  of  reaction  also  increases  due  to  increas¬ 
ing  irreversibilities  as  the  cell  tends  away  from  open  circuit 
conditions.  Fig.  10  shows  the  variation  of  ohmic  and  reactive  heat 
generations  with  cell  voltage.  Once  again,  these  graphs  are  not 
significantly  affected  by  flow  rate  in  the  range  0. 1-1.0  L  min-1 . 
If  there  were  a  more  pronounced  concentration  overpotential 
region  in  the  IV  curve,  then  the  rate  of  increase  of  current  with 
voltage  would  decrease  sharply  at  limiting  conditions,  and  hence 
the  ohmic  heating  contribution  would  be  less  than  those  shown 
in  Fig.  10.  This  would  result  in  lower  temperatures  at  limiting 
conditions,  than  those  predicted  in  Fig.  8.  Nevertheless,  it  is  not 
desirable  to  run  the  fuel  cell  close  to  limiting  conditions.  At 
optimum  conditions,  the  temperature  variations  are  much  more 
manageable. 

4.  Conclusions 

A  three  dimensional  mathematical  model  of  a  PEM  fuel 
cell  equipped  with  a  PBI  membrane  was  developed.  The  model 
accounted  for  all  transport  phenomena  and  polarization  effects. 
Its  predictions  match  well  with  experimental  results  published 
by  Wang  et  al.  [13].  Compared  to  a  previously  developed  two- 
dimensional  model,  the  3D  model  predicted  a  greater  ohmic  con¬ 
tribution  due  to  the  presence  of  ribs,  which  result  in  a  decrease  in 
IV  characteristics.  In  other  words,  the  2D  model  overestimated 
the  performance  by  neglecting  rib  effects. 

Results  show  that  the  greatest  area  of  oxygen  depletion  occur 
in  the  cathode  catalyst  layer  just  under  the  ribs.  This  depletion 
increases  in  the  direction  of  flow,  and  is  more  prominent  at  lower 
supply  gas  flow  rates.  The  temperature  increases  in  the  cell  in  the 
direction  of  the  channel  flow,  as  expected,  since  heat  generated 
is  absorbed  by  the  feed  gases.  The  hottest  points  occur  in  the 
membrane  just  under  the  ribs.  This  is  close  to  the  region  of 
greatest  ohmic  and  reactive  heating.  The  temperature  variation 
across  the  cell  ranges  from  7  to  19  K  at  optimum  conditions. 
At  limiting  conditions,  much  higher  temperature  variations  are 
predicted. 

This  model  did  not  take  into  account  reactant  gas  solubility  at 
the  catalyst  layers,  since  it  assumed  gas  phase  reactions.  Taking 
this  phenomenon  into  account,  as  well  as  dissolved  gas  kinet- 
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ics,  may  render  the  model  better  able  to  predict  mass  limitation 
effects.  This  is  the  subject  of  future  work. 
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